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We propose a unifying rheological framework for dense suspensions of non-Brownian spheres, pre¬ 
dicting the onsets of particle friction and particle inertia as distinct shear thickening mechanisms, 
while capturing quasistatic and soft particle rheology at high volume fractions and shear rates re¬ 
spectively. Discrete element method simulations that take suitable account of hydrodynamic and 
particle-contact interactions corroborate the model predictions, demonstrating both mechanisms of 
shear thickening, and showing that they can occur concurrently with carefully selected particle sur¬ 
face properties under certain flow conditions. Microstructural transitions associated with frictional 
shear thickening are presented. We find very distinctive divergences of both the microstructural and 
dynamic variables with respect to volume fraction in the thickened and non-thickened states. 


I. INTRODUCTION 

Non-Newtonian rheology [l| has been observed and 
studied for centuries in numerous materials, flow regimes 
and ^plications. In this work we focus on shear thicken¬ 
ing [2| in densely packed non-Brownian suspensions of 
bidisperse solid spheres, with and without inertia [3“ 
[^. This rheological phenomenon, in which the shear 
stress required to deform the suspension increases faster 
than linearly with the deformation rate, is regularly 
demonstrated in high volume fraction cornstarch suspen¬ 
sions @ , but is also observed in other particulate systems 
such as dry granular materials at constant volume 
and well characterised model suspensions [l^, and has 
considerable industrial relevance |ll| . The non-Brownian 
limit arises in suspensions of both silica and polymethyl¬ 
methacrylate, for example, under typical shear thicken¬ 
ing conditions [T^ . 

Continuous, linear shear thickening, in which the sus¬ 
pension viscosity is proportional to the shear rate, may 
arise in suspensions below jamming [l3l | when conditions 
are such that particle inertia is relevant [T3 - [T^ . much 
like in dry granular materials m. Other suspensions 
have, however, been observed to shear thicken far more 
severely than in these linear cases, and at Stokes num¬ 
bers considerably less than 1, for which particle inertia 
ought to be negligible [a, flSl - l^ . This behaviour is var¬ 
iously known as “shear jamming”, “dynamic jamming” 
and “discontinuous shear thickening” (DST) [^, and un¬ 
til recently has widely been thought to arise due to either 
the shear-induced formation of “hydroclusters” [l^, [13 > 
mesoscale particle agglomerates stabilised by hydrody¬ 
namic interactions that result in massive dissipation un¬ 
der shear; or dilatancy, the tendency of the suspension to 
increase in volume upon shearing @,[13 and subsequently 
bifurcate into coexisting regimes of inhomogeneous solids 
fraction (l^ - 

A growin g bo dy of experimental [13, [13 compu¬ 
tational [29l - [31j work proyides eyidence that discontinu¬ 
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ous shear thickening can arise because frictional particle- 
particle contacts appear under large loads. Such sus¬ 
pended particles may be either charge stabilised or steri- 
cally stabilised using, for example, polymer hairs grafted 
to the particle surface. Under small loads, the nor¬ 
mal repulsiye forces that arise between particles due to 
this stabilisation are sufficient to preyent direct particle- 
particle contacts, so lubricating layers are maintained. 
Above a critical load P*, the stabilisation is overcome 
and rough particle surfaces come into contact, result¬ 
ing in normal and tangential forces that can be consid¬ 
ered similar to those existing between dry granular par¬ 
ticles [111- The iucreased dissipation resulting from the 
subsequently reorganised microstructure and the tangen¬ 
tial contact forces means very large stresses are required 
to maintain flow. Under this mechanism, the shear thick¬ 
ened state may flow homogeneously, without velocity or 
volume fraction banding |28l| . A rheological model pro¬ 
posed by Wyart and Cates [^ captures this transition 
between frictionless and frictional states, predicting the 
presence of continuous shear thickening at low volume 
fractions, and DST at high volume fractions, where S- 
shaped flow curves could occur with multiple flow states 
existing at a given shear rate but at greatly differing 
stresses. Such flow curves have recently been observed 
experimentally and computationally under im¬ 
posed shear stress. 

In general, the rheology of suspended “nearly-hard” 
spheres can be broadly characterised by the interplay 
between the flow timescale associated with inverse of 
the shear rate I/7, and four competing timescales: a 
Brownian timescale (characterised by the Peclet number 
Pe = inrjf'jcP/ikT, for representative particle diameter 
d and interstitial fluid viscosity rjf); a timescale associ¬ 
ated with the stabilising repulsion (numerically this has 
been referred to as I/70 = ^'Kr]fd?/F'^^, for repulsive 
force magnitude [sO]; an inertial timescale (char¬ 

acterised by the Stokes number St = pd'^j/rjf, where p 
is a representative suspension density and St = 1 de¬ 
lineates viscous and inertial flows) and a timescale as¬ 
sociated with the stiffness of the particles (for example 
d/ \/kn/ pd, where kn is related to the Young’s modulus 
of the particles) Q • So far, these diverse scales have not 
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FIG. 1. Steady state rheological regime map for a shear 
thickening suspension, illustrating the frictional thickening 
transition within the viscous regime, shear jamming, inertial 
shear thickening, quasistatic behaviour and deformational be¬ 
haviour associated with soft particle rheology. 

been probed simultaneously in a single suspension either 
experimentally or computationally, though there are nu¬ 
merous recent examples of transitions across regimes that 
suggest they represent isolated re gion s of a single rheo¬ 
logical map for dense suspensions |l3 . 1 ^ . 1 ^ . . 

In the present work, we focus on the non-Brownian 
limit (i.e. Pe —>■ oo), and demonstrate that constitutive 
models proposed for shear thickening in the non-inertial 
limit |33l| and for capturing the non-inertial (viscous) to 
inertial transition |l5l Il6l| . can be unified to place fric¬ 
tional shear thickening in the wider context of dense sus¬ 
pension rheological regimes, also accounting for the ef¬ 
fects of finite particle hardness. We then perform discrete 
element method 1391 l40ll simulations combining hydrody¬ 
namic lubrication |4l| with a suitable particle-particle 
interaction model proposed by Mari et al. that 

can capture the bulk steady-state rheological behaviour 
associated with frictional shear thickening under imposed 
shear rate. We demonstrate that the timescale for fric¬ 
tional shear thickening can be made to coincide with that 
for inertial shear thickening by careful tuning of parti¬ 
cle surface properties, hinting at novel suspension flow 
curves that have yet to be observed experimentally. Fi¬ 
nally, we highlight microstructural properties associated 
with the frictional thickening transition, identifying very 
well defined structural and dynamic signatures that may 
prove useful in interpretation and analysis of future rheo- 
imaging data for shear thickening suspensions. 

II. CONSTITUTIVE MODELLING AND FLOW 
REGIME MAP 

We hrst present a rheological equation that is able 
to capture the viscous, inertial, quasistatic and soft- 


particle flow regimes [l^ . This model is inspired by 
the inertial number model i i (for inertial number 
I I = 'jdl y/Pjp, with confining pressure P), its extension 
to viscous flows [i^ (for viscous number ly = ivt/P), 
and their proposed unification by Trulsson [15l |. The 
equation gives a prediction for the scaled (by particle 
hardness) pressure (P = Pd/kn) as a function of the 
scaled shear rate (7 = jd/y/kn/pd) and the departure 
of the solids volume fraction 0 from its critical value 
for jamming [l^ . (f>c, in each of three regimes: 1 ) the 
hard particle regime corresponding to viscous and iner¬ 
tial flows; 2 ) the soft particle regime corresponding to 
deformable particle flows; 3) the quasistatic, “jammed” 
regime: 


ahardl^-</'c| +aLrd?7/l^-</'c| 

^ 7 , (la) 

contact fluid 


Psoft ^softT ^soft^/T ; 

contact fluid 

(lb) 

PqS = OiQslp - , 

contact 

(Ic) 


with the constants given by Ness and Sun [T^. An 
arbitrarv blending function is chosen, following Chialvo 
et al. [7|, to combine pressure predictions from each of 
the expressions. The corresponding shear stresses a^y 
are obtained from a p{K = Iv + a/|) model [SSl’ 
extension of the commonplace //(//) rheology [9|: 

^{K) = Pi + 1.2A:^/2 -b 0.5Ab (2) 

The rheology predicted by Equations [T] and [2] is pre¬ 
sented in detail in Ref. [1^. We next incorporate a fric¬ 
tional shear thickening mechanism into this model using 
a stress-dependent critical volume fraction (t>r, following 
the approach used by Wyart and Cates [^. A stress- 
dependent effective friction is introduced, recognising 
that the critical volume fraction depends on the interpar¬ 
ticle friction coefficient Pp [i^, varying from </>„ ~ 0.58 
to (po Ri 0.64 in the limits of frictional {pp = 1) and fric¬ 
tionless {pp = 0) particles, respectively. From the simu¬ 
lation model described in Sec IIII A[ we find that under 
shear flow, the rescaled pairwise particle-particle contact 
force magnitudes 0 = |FF|/Pd^ are distributed accord¬ 
ing to 

PDF(0) = a(l — 6 exp(—0^)) exp(—c0), (3) 

consistent with previous authors [3 The fraction of 
particle contacts for which the repulsive force magnitude 
is exceed and friction is activated is therefore given 
by 

POO 

/= / PDF(0)d0, (4) 

JFCL/pd2 
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implying (except for very weak contacts) that frictional 
forces arise in the system above P* according to / oc 
exp(—P*/P). We therefore use / to represent a transi¬ 
tion from frictionless to frictional rheology with increas¬ 
ing P. The value set for P* (which is directly related 
to the repulsive force magnitude F'^^) determines the 
critical pressure (or critical characteristic shear rate) at 
which the model will begin to predict frictional rheology, 
as described later. We subsequently calculate the (stress- 
dependent) critical volume fraction for jamming (pc using 
an expression similar to the crossover function proposed 
by Wyart and Cates 

(pc = (pmf + (po{^ - !)■ ( 5 ) 

The expression for along with that proposed 

by dll, assumes a constant value for the macroscopic 
friction fii (= <Jxy/P = 0.38) in the limit of —>■ 0. 
As demonstrated by da Cruz [^, fii is actually strongly 
dependent on interparticle friction /ip, particularly for 
fip close to 0, so assuming a constant value across shear 
thickening states is clearly not appropriate. We therefore 
propose a similar crossover function for /ii (consistent 
with that proposed by Sun and Sundaresan for dry 
granular materials) which we find gives excellent agree¬ 
ment with the following simulation results 

+ -/)> ( 6 ) 

where /ii^ = 0.41 and /iig =0.11 [d^. We obtain the vis¬ 
cosity of the suspension relative to that of the interstitial 
fluid according to r]s = 

We obtain the flow curves presented in Fig [1] from 
Eqs. [U [21 El and El Below pc, the model predicts vis¬ 
cous rheology for Stokes numbers less than unity, iner¬ 
tial flow at higher Stokes numbers, and “intermediate” 
rheology at extremely high shear rates or for soft parti¬ 
cle suspensions (i.e. emulsions) where large deformations 
are possible (strictly, “soft” particle rheology is ob¬ 
served when 7 —?► dy/kn/pd) Q- Above pc, quasistatic 
rheology is observed for low and moderate Stokes num¬ 
bers, with a tendency towards soft particle rheology at 
very high rates. The addition of pressure dependence in 
pc gives rise to the frictional thickening and S-shaped 
flow behaviour, as shown in Fig [T1 within the viscous 
flow regime, and the hypothetical shear jamming transi¬ 
tion that may occur between the viscous and quasistatic 
regimes for particles of finite stiffness. The fact that the 
onset stress P* varies with particle contact properties 
(specifically ^ the repulsive force magnitude) implies 
that the transition to frictional behaviour might occur at 
different regions of the flow map. We demonstrate this 
behaviour in Fig|2l by increasing (from Fig|2](a) to (c)) 
the magnitude of the onset stress P* in the definition of 
/, delaying the onset of the frictional behaviour governed 
by Eqs. Eland El such that it occurs at higher Stokes num¬ 
bers. We obtain the same characteristic frictional shear 
thickening flow curve predictions for each of Figs E^-c, 
with the added phenomena of inertia appearing at a pre¬ 
scribed point in relation to P* (inEb andEt)- Such shear 


thickening will be demonstrated by particle simulations 
in the next section. 


III. SHEAR FLOW SIMULATIONS 


A. Numerical method 


The equations of motion for non-Brownian particles 
suspended in a fluid can be written simply as [47( | 


4(:)=E(r). (^) 


for particles of mass m with translational and rotational 
velocity vectors v and w respectively, subjected to force 
and torque vectors F and F respectively. In this work we 
limit the forces and torques to those arising due to di¬ 
rect particle contacts (F“, F“) and those arising through 
hydrodynamic interactions (F^,F^). Full solution of the 
pairwise hydrodynamic forces has traditionally been done 
using the Stokesian Dynamics algorithm [^, [^ , though 
its great computational expense makes large (or very 
dense) simulations challenging. For highly packed sus¬ 
pensions, the divergent lubrication resistances that arise 
between extremely close particles dominate the hydrody¬ 
namic interaction, so F^, F^ can be approximated by sum¬ 
ming pairwise lubrication forces among nearest neigh¬ 
bouring particles [3 E Ham 113 ■ For an interaction 
between particles i and j, the force and torque due to 
hydrodynamics can therefore be expressed as 


F^- - 

= -asq6TTr]f{v^ - 

w) •' 


fScl) 


- ash&T^rjfi^i - Vj 

)•(!- 

- riijiiij), 


r^.. : 

= -apuTTrjfdf^LJ^ 


■ (I - nyriy) 





d- 

(8b) 



— 

Y (W, X Fy . 



where is the vector pointing from particle j to parti¬ 
cle i, and with squeeze shear ash and p ump apu resis¬ 
tance terms as derived by Kim and Karrila[5l| and given 
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in Eq[9]for particle diameters di and dj, with /3 = dj/df. 


dl l + 7/3 + /3^ d. 

(1 +/3)2 2/leff 5(l + ^)3 2 ""Vs/iesy 


1 + 18/3 - 29/3^ + 18/33 +/34 ^2 


21(l + /3)4 


^ Q 2 + /3 + 2/3^ di ( d, 
ash = 4d — in 


In 


4/ieff v2/leff 


di_ \ 

ef 


(9a) 


+ 4 


15(l + /3)3 2“‘V2deff, 

16 - 45/3 + 58/3^ - 45/33 + 16/3^ df 


375(1 + /3)4 


npit — /3“ 


4 + /3 


■In 


d,; 


^10(l+/3)2‘“V2deff, 

32 - 33/3 + 83/32+43/33 d. 


250/33 


■ In 


4deff 


/ dj 


In 


2d0ff V2d0ff 


dj 

2 d 0 ff y 

(9b) 


(9c) 


For each pairwise interaction, the surface-to-surface dis¬ 
tance, d, is calculated according to d = |rij| — for 

centre-to-centre vector Tij. Recent experimental |^ . l3^ 
and computational [12, [S^ work indicates that direct 
particle-particle contacts play a significant role in de¬ 
termining steady state paste viscosity. To permit such 
contacts in the present model, we truncate the lubri¬ 
cation divergence and regularize the contact singularity 
at a typical asperity length scale dmin = O.OOldij (for 
weighted average particle diameter d^ = ^ ), i.e., set¬ 

ting h = hjnin in the force calculation, when h < 

The effective interparticle gap used in the force calcula¬ 
tion, hes, is therefore given by 


defF — 


for h > d-min 
min otherwise. 


( 10 ) 


For computational efficiency, the lubrication forces are 
omitted when the interparticle gap d is greater than 
dmax = 0.05dij. The volume fraction is sufficiently high 
that all particles have numerous neighbours within this 
range, so such an omission is inconsequential to the dy¬ 
namics. When the lubrication force is overcome and par¬ 
ticle surfaces come into contact (this occurs when h < 0), 
their interaction is defined according to a linear spring 
model [1^, with normal (F“’” ) and tangential (F“’‘ ) 
force and torque F"^ given by 




(11a) 


F'rf = -fctun 


1J5 


(11b) 


n = -|(ny X F^‘), (11c) 

for a collision between particles i and j with normal and 
tangential spring stiffnesses and kf respectively, parti¬ 
cle overlap S (equal to —d) and tangential displacement 


Uij. We note that the damping arising from the hydro¬ 
dynamics is always sufficient to achieve a steady state 
without employing a thermostat, and further damping 
in the particle-contact model is omitted for simplicity. 

We employ the Critical Load Model (CLM) for inter¬ 
particle friction [2^ [s^l, to distinguish between weakly 
interacting particles, those that interact via the normal 
force deriving from stabilisation, and strongly interacting 
particles, those whose surfaces come into frictional con¬ 
tact. This model gives an additional stress scale for the 
particle interaction, which, numerically, is the origin of 
the onset stress for shear thickening P*. An interparti¬ 
cle Coulomb friction coefficient is defined according to 
iFi’jl < setting a maximum value for the tan¬ 

gential force exerted during a collision. For each pairwise 
collision, the value of fip is dependent upon the normal 
force between the interacting particles and some critical 
normal force magnitude representing the magni¬ 

tude of the stabilising repulsive force, such that 




1 for \F‘r’J\ > 
0 otherwise 


( 12 ) 


As a result of the CLM, particles that interact through 
weak forces, i.e. collisions where (5—5-0, are friction¬ 
less, while interactions with large normal forces are fric¬ 
tional. This particle potential represents a physical sce¬ 
nario closer to electrostatic rather than polymer hair 
driven normal repulsion. The particle overlaps required 
to exceed are, at their absolute maximum, of order 
lO-^dij. 

Isotropic particle assemblies are generated in a 3- 
dimensional periodic domain of volume V. It is deter¬ 
mined that approximately 5000 spheres are sufficient to 
capture the bulk rheology and microstructural phenom¬ 
ena independently of the system size. Bidispersity at a 
diameter ratio of I : 1.4 and volume ratio of 1 : 1 is used 
to minimize crystallization during flow [3^ . The parti¬ 
cle assembly is sheared to steady state at constant rate 
7 and constant volume, equivalent to the application of 
Lees-Fdwards boundary conditions [s^. The bulk stress, 
decomposed into contributions due to the hydrodynamic 
interaction and the particle-particle interaction, is cal¬ 
culated from the particle force data [l^, and given by 
Fqs. I13al and Il3bl 


1 

(I3a) 


(13b) 


i 


Data from 20 realizations with randomized initial particle 
positions are used to obtain ensemble-averaged stresses, 
which are further averaged over time in the steady-state. 
Under simple shear flow, the relevant stresses that will 
be discussed are the xy component axy (= '^xy'^^xy)^ 
flow direction x and gradient direction y, and the mean 
normal stress P = \{(Jxx+(^yy + (^zz), for axx = cTxx + ^^xx 
etc. 
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FIG. 2. Shear thickening transition for three values of P* 
(or, equivalently, (a) Frictional shear thickening occurs 

in the absence of inertia. Dashed line and arrow demonstrates 
the relative location of rheological data presented by Ness and 
Sun [ 13 ] : (b) Frictional shear thickening occurs concurrently 
with the onset of inertia; (c) Frictional shear thickening occurs 
in the presence of inertia. Coloured circles represent discrete 
element method simulation results; solid black lines represent 
constitutive model predictions; dotted black lines represent 
P*. 


B. Macroscopic flow behaviour 

Transitions between the viscous, inertial, soft par¬ 
ticle and quasistatic regimes, as they are depicted in 
Figlll have been previously captured by discrete element 
method simulations and well characterised [l^. The 
steady state shear thickening behaviour predicted by the 
simulation model described in Sec IIII Al is presented as 
solid coloured symbols in Fig[2l We first focus on the re¬ 
sults in Fig[ 2 ^, which correspond directly to the frictional 
thickening transition highlighted within the viscous flow 
regime in Fig[TJ Following [^, the shear rate 7 is scaled 
with the reciprocal of a characteristic timescale for the 
relaxation of a frictional contact in a viscous fluid, given 
by 7 o = /^nrifd?. Consistent with the results ob¬ 
tained by Mari et al. [ 2 ^ [s^ , shear thickening between 
two quasi-Newtonian flow regimes is observed to occur 
at an onset stress P*, independent of volume fraction 
and given by the dashed black line in Fig [2^. Far below 
the onset stress, particles interact through forces pre¬ 
dominantly |F^’"| < that is, the forces are not 

sufficiently large to overcome the stabilisation, so fric¬ 
tional particle surfaces do not often come into contact. 
Conversely, the stabilisation is nearly always overcome 
(so contacts are nearly always frictional) at stresses far 
above the onset stress. We therefore make a distinction 
between purely frictionless behaviour at 7/70 = 0.01 and 
purely frictional behaviour at 7/70 = 1. The solid black 
lines represent predictions from the constitutive model 
described previously. The value of the onset stress is de¬ 
termined by the magnitude of specified in the con¬ 
tact potential, and is inferred from the simulation data. 
The annotation in Fig [5^ illustrates the relative position 
of the simulation data presented by Ness and Sun [l3| . 

For volume fractions below approximately (f) = 0.53, 
the rheology exhibits continuous shear thickening be¬ 
haviour, while between RiO.53 and ~0.58 the thickening is 
discontinuous, in that the constitutive model flow curves 
(solid black lines) exhibit the characteristic ‘S-shaped’ 
phenomena. The simulation data do not populate the 
‘S-shaped’ region, probably because the simulations were 
performed for steady states at enforced constant shear 
rate, while the nature of flow in this regime is high un¬ 
stable in reality. Probing the ‘S-shaped’ region through 
other simulation protocols is the subject of ongoing inves¬ 
tigation and will be reported on in the future. For volume 
fractions above 4>m, the material “shear-jams” above P*, 
as illustrated in Fig [T] When the onset stress is exceeded 
above (j)m, the material transitions from below to above 
its critical volume fraction, meaning the flow moves from 
a flowing, viscous state to a jammed state. Experimen¬ 
tally, this may be manifested as complete flow cessation, 
surface fracture, microstructural inhomogeneity, or vol¬ 
ume fraction bifurcation, depending on the nature of the 
rheometer. Particle overlaps are allowed in the simula¬ 
tions, so the flow can enter a quasistatic state above jam¬ 
ming, in which the shear stress is rate-independent [l6|. 

To bring the frictional thickening transition nearer 
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to the inertial regime computationally, we simply re¬ 
duce the viscosity of the interstitial fluid rjf, modifying 
7o(= /^TT'qfdlj) and effectively moving the 0.01 < 
7/70 < 0-1 window to a higher range of Stokes numbers. 
We can achieve an analogous effect by adjusting , 
comparable to modifying either the particle size or the 
particle surface chemistry experimentally. In terms of 
shear thickening, the effect of this adjustment is to alter 
the magnitude of the onset stress for frictional contacts 
such that it occurs in the vicinity of any desired Stokes 
number. Flow curves are presented for an onset stress 
that occurs close to St= 1, Fig[2b, and for an onset stress 
that occurs at very high St, Fig[2j:. In each case, a tran¬ 
sition is observed between the frictionless and frictional 
states, similarly to the totally viscous (St <C 0) case. 
In Fig the frictionless regime is observed for Stokes 
numbers up to around unity. Below this point, the sus¬ 
pension viscosity is independent of the Stokes number. 
For larger Stokes numbers, we observe frictional, inertial 
flow, with (Jxyldfi oc 7. The linear scaling of viscosity 
with shear rate above St= 1 is due to inertial effects; the 
larger jump in viscosity (i.e. the super-linear behaviour 
between 7/70 = 0.1 and 7/70 = 1) is due to the onset 
of frictional contacts. We verify that the flow in each 
of these limits remains frictionless and frictional respec¬ 
tively at the microscale by examining the fraction of par¬ 
ticle contacts that have exceeded F^^. It is noted that 
the result in Fig [ 2)3 corresponds directly to the shear 
thickening phenomenon observed in the simulations re¬ 
ported by Fernandez [^, with a low shear rate regime 
in which lubrication dominates and particle friction is 
absent and a high shear regime dominated by friction 
with a viscosity that depends linearly on the shear rate. 
Fernandez also reports experimental findings for shear 
flow of polymer coated quartz miroparticle suspensions 
that appear qualitatively similar to Fig [ 2 ) 3 , though it is 
not clear whether the Stokes number is appropriate for 
such a comparison. Indeed, a similar set of experimen¬ 
tal findings [sj were previously attributed to enhanced 
dynamic friction due to increased resistance to fluid flow 
in the polymer layer. In Fig [Jt, the onset stress occurs 
at St ss 300, so both the frictionless and frictional states 
(again, these are verified by appealing to the frictional 
of individual contacts) exhibit o'xy/rjf'y oc 7 scaling, with 
super-linear behaviour representing the frictional transi¬ 
tion. In each case, we obtain from the simulations excel¬ 
lent agreement with theoretical predictions in the limits 
of fully frictionless and fully frictional flow. 

These novel flow curves clearly illustrate the distinc¬ 
tion between shear thickening driven by friction and by 
inertia. By tuning the particle and fluid properties ap¬ 
propriately, we have demonstrated that although (ex¬ 
perimentally) the regimes may seem highly distinct and 
therefore challenging to achieve within a single system, 
both mechanisms might be made to arise concurrently, 
giving rise to new rheological behaviour. The challenge 
remains to achieve a sufficient understanding of the roles 
of particle surface chemistry, particle size and suspending 



<!> 

(a) 



4 ^ 

(b) 


FIG. 3 . (a) Divergence of viscosity contributions and model 

prediction; (b) Divergences of normal and tangential contact 
forces. 


fluid properties to realise and utilise these flow regimes 
experimentally. 

For the entirely viscous case presented in Fig [2ja,, we 
isolate the contact and fluid contributions to the viscos¬ 
ity and plot them against volume fraction for the fric¬ 
tionless (7/70 = 0.01) and frictional (7/70 = 1) limits in 
Fig [3^. It is noted that analogous results are obtained 
for the inertia-dominated cases, though the magnitude 
of the viscosity is increased consistent with the linear 
viscosity scalings demonstrated in Figs[2l3-c. Comparing 
the jumps from the frictionless to the frictional branch, it 
is demonstrated that the main increase in viscosity upon 
shear thickening is due to the contact contribution, while 
there is only a very minor increase in the fluid contribu¬ 
tion (appreciable only at high volume fractions). While 
this suggests a configurational change leading to a change 
in the mean fluid film thickness or film number (result¬ 
ing in a slightly increased fluid stress), it is not consis¬ 
tent with the notion of a large macroscopic transition 
to hydroclustering and a corresponding massive in¬ 
crease in viscous dissipation. We further decompose the 
contact stress into normal and tangential components, 
Fig [ 3 ) 3 . We find that although the major difference be¬ 
tween the frictionless and frictional limits at the indi¬ 
vidual particle level is the presence of tangential contact 
forces, the main contributor to the increase in the con¬ 
tact stress is in fact the normal component, rather than 
the tangential component, further corroborating the ma¬ 
jor role played by the particle configuration change in- 
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duced by friction. This change can be understood from 
the perspective that the available degrees of freedom for 
particle motion decrease at the onset of frictional con¬ 
tacts, in that frictional particle assemblies require four 
contacts per particle for mechanical stability, while fric¬ 
tionless ones require six [s^. At fixed volume fraction, 
the transition to frictional behaviour is therefore mani¬ 
fested as an increased resistance to flow that necessitates 
greater particle overlaps and results in higher particle 
pressure. Interestingly, reducing the available degrees 
of freedom by means other than particle friction leads 
to the same observation. In a separate simulation we 
model steady shear at 7/70 = 0.01, and at time ti we 
set all z—components of particle velocity and forces to 
zero, effectively imposing a 2D flow constraint. A large 
increase in the particle contribution to the stress is ob¬ 
served, consistent with the shear thickening behaviour 
presented here, that dissipates at a later time t 2 when 
the 2D flow constraint is relaxed. The dominant role 
of contacts, and the sensitivity to their nature (whether 
tangential forces may be sustained in addition to normal 
forces), remains a contentious issue; these results add fur¬ 
ther weight to the argument for frictional contacts as a 
crucial contributor to (non-inertial) shear thickening. 

We have therefore demonstrated that the frictional 
thickening mechanism, modelled via a load-dependent 
particle friction in DEM and captured by the constitutive 
model, can occur within a variety of flow regimes and may 
thus couple or compete with inertial thickening. Both the 
DEM simulation and the constitutive model capture con¬ 
sistently such shear thickening phenomena. By isolating 
the contact and hydrodynamic contributions to the shear 
stress, we have shown that the shear thickening transi¬ 
tion is heavily dominated by particle contacts as opposed 
to hydrodynamic effects. The large increase of contact 
stress upon shear thickening has been attributed mainly 
to normal contact forces, though we note that the tangen¬ 
tial contact forces present in the thickened regime largely 
exceed the normal forces present in the non-thickened 
state. 


C. Microscopic analysis 

We use microscale information to further characterise 
the distinction between the thickened (frictional) and 
non-thickened (frictionless) states. The microstructure 
is quantified using the particle-particle contact number 
and the fabric, or net-orientation, of these contacts, while 
particle-level dynamics are quantified using correlation 
functions in displacement and velocity. Work in rheo- 
imaging of colloidal systems [H, [s^l has demonstrated 
the potential for these dynamic properties to be obtained 
and quantified experimentally for shear thickening mate¬ 
rials. In addition, further experiments have led to 
their successful quantification for a sheared, highly poly- 
disperse emulsion, using confocal imaging. We anticipate 
that future such analyses for shear thickening suspensions 


will benefit from, and corroborate, the results presented 
in this paper. 

The microstructure is characterised based on two sep¬ 
arate length scales. For the contact number Z, we calcu¬ 
late (a) the average number of frictionless, 
interactions per particle; (b) the average number of fric¬ 
tional, |F°’"| > contacts per particle. For the con¬ 
tact fabric, we adopt the formulation used by Sun and 
Sundaresan [i^ l 



Nc 



(14) 


Under shear flow, contacts preferentially align along 
the compressive axis at (or close to) 45° to the x— 
and y—axes (the flow and gradient directions respec¬ 
tively), with the corresponding shear component of A, 
IA 12 I quantifying the extent of the anisotropic alignment. 
Ai 2 = 0.5 represents perfect alignment of all contacts in 
the compressive axis, while A 12 = 0 represents perfect 
isotropy. As with the contact number, we quantify the 
shear component of the fabric based on both the network 
of frictionless contacts and that of frictional contacts. 
The non-afSne motion is quantified by first obtaining a 
coarse-grained velocity profile for the shearing flow and 
subtracting the appropriate value of this coarse-grained 
velocity from each particle’s velocity vector to obtain the 
“fluctuating” velocity component. From these fluctuat¬ 
ing velocities we obtain the mean squared displacement 
(MSD), averaged across particles and time steps. Plot¬ 
ting the MSD ((cc^)) versus strain magnitude (jt) yields 
a linear diffusive behaviour that is independent of Stokes 
number, for the case of zero inertia. We take the gradi¬ 
ent of this slope D^, as an effective diffusion coefficient 
{i.e. = Dxjt), though strictly the diffusion coeffi¬ 

cient in this case is T>xi- In addition, we calculate the 
correlation of the fluctuating velocity vectors according 
to Lois [illillil 


Ei EiM'^dUjI -r) 


(15) 


where W, Vj are particle velocity vectors averaged over a 
length of time sufficient to give an averaged particle dis¬ 
placement due to the mean flow of approximately O.hd. 
It is found that the correlation decays approximately ex¬ 
ponentially with the distance between particle centres r. 
We therefore fit a simple functional form C'(r) = 
where ^ takes units of particle diameter and is hereafter 
referred to as the “correlation length”, and fc is a constant 
prefactor. The evolution of these microscale quantities is 
presented in Fig |4] as a function of 7/70 at three vol¬ 
ume fractions, and as a function of volume fraction at 
7/70 = 0-01 (frictionless) and 7/70 = 1 (frictional), for 
the case of zero inertia. 

In Figd^ we demonstrate the increasing number of fric¬ 
tional contacts and the diminishing of frictionless inter¬ 
actions as the shear rate is increased and the onset stress 





FIG. 4 . (a)-(d) Evolution of microscale structures and dynamics across the frictional shear-thickening transition, (a) Mean 

number of particle-particle contacts. Squares represent low force contacts for which friction is not activated; circles represent 
high force contacts for which friction is activated, (b) Shear component of the fabric tensor. Squares represent low force contacts 
for which friction is not activated; circles represent high force contacts for which friction is activated, (c) Effective diffusion 
coefficient, (d) Velocity correlation length as defined by Eg 1151 (e)-(h) Evolution of microscale structures and dynamics with 
volume fraction, for frictional and non-frictional states, (e) Mean number of particle-particle contacts. Stars represent low 
force contacts for which friction is not activated; circles represent high force contacts for which friction is activated, (f) Shear 
component of the fabric tensor. Stars represent low force contacts for which friction is not activated; circles represent high 
force contacts for which friction is activated, (g) Effective diffusion coefficient, (h) Velocity correlation length. 


is exceeded. The results presented here are consistent 
with the evolution of the fraction of frictional contacts 
presented by Mari et al. [s^. It is further noted that for 
a fixed volume fraction, the number of direct particle- 
particle contacts that exist in the frictional, shear thick¬ 
ened state exceeds the number of frictionless, normal in¬ 
teractions that were present in the non-thickened state. 
This suggests that in the process of becoming frictional, 
the particles have rearranged into a distinct microstruc- 
tural configuration. The evolution of A 12 , Fig[4)D, con¬ 
firms this. In the non-thickened state, we observe dis¬ 
tinct microstructures for the networks of frictionless and 
frictional contacts, though the number of frictional con¬ 
tacts is very small. At each volume fraction, contacts for 
which is exceeded tend to be aligned more strongly 
with the compressive flow direction than those contacts 
for which friction is not activated. Upon shear thicken¬ 
ing, however, the fabric of the frictional contact network 
moves closer to zero, while the frictionless fabric disap¬ 
pears due to an absence of such interactions far above 


the onset stress. The microstructural information sug¬ 
gests that the shear thickening transition brings the par¬ 
ticle configuration closer to what might be expected for 
a quasistatic, rate-independent flow, where Z Ri 4 and 
Ai 2 « —0.03 [ 4 ^ [ 5 ^ for the frictional cases. Shear thick¬ 
ening can therefore, in this sense, be considered analo¬ 
gous to an increase in volume fraction at constant fric¬ 
tion, in that the central change in each case is that the 
departure from the critical volume fraction, quantified as 
\(j) — decreases. It is noted that the inertial cases 
exhibit very similar behaviour with respect to the mi¬ 
crostructural properties. The exception is a very modest 
increase in contact number, smaller than 10%, for the 
inertia-dominated flows. 

Particle level dynamics are found to exhibit analogous 
behaviour across the thickening transition. In Fig 0);, 
we show that for each volume fraction studied, there is 
a significant jump in the effective diffusion as the sus¬ 
pension shear thickens. As with the microstructure, this 
is consistent with the suspension becoming closer to its 


























9 


jamming volume fraction as the shear rate or stress is in¬ 
creased. As the extent of frictional behaviour increases, 
particles form more contacts and are required to deviate 
further from an affine trajectory in order to pass each 
other as they are subjected to shear at the same or higher 
rate. The average displacement deviation from the mean 
flow therefore undergoes a step change. A similar transi¬ 
tion is observed in the velocity correlation length, FigHJi, 
which indicates that in the shear thickened regime, par¬ 
ticle trajectories tend to be more correlated with those 
of their immediate neighbours, suggesting a tendency to¬ 
wards collective motion of particle groups. Though this 
is, indeed, qualitatively reminiscent is some respects to 
the “hydroclustering” behaviour previously reported [13 j 
we note that the particle groups that collectively move 
in the present simulations are found to be unanimously 
under frictional contact, rather than separated by lubri¬ 
cation layers. We therefore suggest that while clustering 
is apparent, it is, in this case, more accurately described 
as frictional- rather than hydro-clustering. Indeed, ex¬ 
perimental techniques that purportedly demonstrate hy¬ 
droclustering are in fact unable to distinguish between 
these two mechanisms of dynamic correlation [^ . 


We demonstrate in Figldjo that each of the transitions 
presented in Fig [4^ corresponds to a shift between dis¬ 
tinct and well-defined branches in each of the microscale 
parameters investigated, in the thickened ( 7/70 = 1 ) and 
non-thickened ( 7/70 = 0.01) limits. In addition, it is 
observed that the differing divergences of each pair of 
branches with volume fraction is consistent with that 
observed for the bulk suspension viscosity, namely the 
frictional, high stress branch diverges near (/) = (fm and 
is thereafter consistent with quasistatic behavior (d^ . 
while the frictionless, low stress branch diverges towards 
4> = (t>RCP- We note that the divergence is less clear 
for the correlation length f than for the other microscale 
parameters investigated, which can be attributed to the 
relatively small domain size, which places limits on the 
length scale over which correlations can be observed. The 
distinction between branches, however, is convincing. A 
demonstration of such contrasting microscale divergences 
in the thickened and non-thickened rheological regimes in 
an experimental system would prove invaluable in corrob¬ 
orating this work, and in highlighting the essential role 
of friction as the origin of the distinct rheologies below 


and above shear thickening. 

IV. CONCLUDING REMARKS 

In this paper we have explicitly demonstrated the dis¬ 
tinction between frictional and inertial shear thickening 
mechanisms, and illustrated their presence as separate 
regimes on the broader flow map of dense suspensions. In 
practice, frictional shear thickening is typically observed 
for colloidal (d ^ I ^m) suspensions for which inertia 
is always absent (or negligible), while the inertial shear 
thickening typically occurs in granular d > 100 /rm) sus¬ 
pensions for which friction is always present (or more 
accurately starting from exceedingly low Stokes num¬ 
bers). Our simulation results suggest that in principle 
thickening may occur in a mixed mode with both mech¬ 
anisms playing a role. This may indeed be the case for 
a suspension of mixed colloidal and granular particles, 
as hinted by the recent experiments on shear thickening 
with intermediate particle sizes [s^ . which indicate that 
the frictional thickening onset stress scales with the in¬ 
verse square of particle size. There are of course many 
other possible scenarios where mixed thickening could oc¬ 
cur as suggested from our simulations, though it remains 
to be seen whether such rheology can be observed exper¬ 
imentally. 

Transitions in the microstructural and dynamic vari¬ 
ables are observed across the frictional thickening transi¬ 
tion, and we have shown that the microstructure of the 
shear thickened and non-thickened states exhibit distinct 
divergences with respect to volume fraction, indicating 
the microscale mechanism for the same behaviour of the 
bulk suspension viscosity. We expect the results pre¬ 
sented here to provide useful means of analysing new re¬ 
sults obtained from particle microscopy and imaging of 
model experimental systems. 
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